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Abstract 

The solution of a (stochastic) differential equation can be locally approxi- 
mated by a (stochastic) expansion. If the vector field of the differential equation 
is a polynomial, the corresponding expansion is a linear combination of iterated 
integrals of the drivers and can be calculated using Picard Iterations. However, 
such expansions grow exponentially fast in their number of terms, due to their 
specific algebra, rendering their practical use limited. 

We present a Mathematica procedure that addresses this issue by re-parametrising 
the polynomials and distributing the load in as small as possible parts that can 
be processed and manipulated independently, thus alleviating large memory 
requirements and being perfectly suited for parallelized computation. We also 
present an iterative implementation of the shuffle product (as opposed to a 
recursive one, more usually implemented) as well as a fast way for calculating 
the expectation of iterated Stratonovich integrals for Brownian Motion. 
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1 Motivation and mathematical background 

In this section, we introduce the mathematical background and motivation for 
manipulating expansions and iterated integrals. The next subsection introduces 
the Picard procedure, a simple iterative way to derive local approximation of 
the solution of a differential equation. Iterated integrals are then introduced 
and the two are finally combined to define the expansions. 
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1.1 Motivation and notation 

Consider the following differential equation: 

dY t = f(Y t ,0)dX t (1) 

where Y t £ M m , X t £ R n and f(.,6) G R mxn . The parameters of the 
function / are collected in the variable 9. We call each j'*' a driver of the 
differential equation. We assume the functions fji, i E {1 . . . n}, j 6 {1... m} 
to be polynomials. The initial value of Yt is set to Yq. 

The objective is to derive a local approximation of the solution in terms of 
/, Yq and the iterated integrals of Yt, i.e. that is integrals of the form 

dxi\...dxn 

0<ui<---<n fe <T 

for any word r = (ri,...,Tfc) constituted of the letter n 6 {l,...,n}, i = 
1, . . . j k.. Such expansions play an important role in the theory of rough paths, 
allowing one to define such a differential equation for a large class of drivers X 
( [3] ) . They have also been used recently for parameter estimation of SDE ( [4 J ) . 

1.2 Picard iterations 

Picard iterations provide a way for deriving local approximations of solutions 
of differential equation. They are defined as: 

Y , T (0) = Y -Y = (2) 

y ,T(r + l,i) = Yl J fiA Y '(r,j))dX® (3) 



i=l 



(4) 



where Y^t^j) is the j th component of the approximation of lo,r = Jq dY s 
Yt — Yq after r iterations. Thus, the first iteration gives: 



Y , T (l,j)= f fj,i(Y )dxP + ...+ I f jtn (Y )dXV) 
Jo Jo 



(5) 



Note that if we are interested in the actual value of the solution at a time 
t € [0,T], its Picard approximation is given by 

Y t (r,j)=Y + Y 0tt (r,j). 

One of the main successes in the theory of rough paths was to give precise 
conditions on X and / for Picard iterations to converge (|3j). 

1.3 Iterated integrals 

Iterated integrals are integrals of the form: 

*#= /•••/ dX^...dXif 

J J S<U\<...<Uk<t 
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where r = (n, . . . , t/~) is called a word, with letters t% 6 {1 . . . n}, i = 1, . . . , k. 
By definition, integrating an iterated integral produces an iterated integral: 



r x^ s dx^ = r [...[ dx^...dx^dx(j) ( 6 ) 

JO JO J Jo<ui<...<u k <s 

= xfe*-^ (7) 

If X is a geometric p-rough path, i.e. it can be approximated by paths of 
bounded variation (see [3] for a precise definition), then the integrals obey the 
usual integration-by-parts rule, which can be generalized as follows: 

QSrUp 

The shuffle product U of two words r and p is the set of all words using the 
letters in r and p such that they are in their original order. For example 
13U42 = {1342,1432,4132,1423,4213,4132} but 1324, for example, does not 
belong to the set. 

This result holds only in the case of deterministic drivers X, or Stratonovitch 
integrals. A similar relation exists for Ito integrals but requires a small cor- 
recting term ([5]) making them slightly less practical in this context. A simple 
transformation allows diffusion processes to be defined with respect to ltd or 
Stratonovitch integrals (PQ): 

dyt = p{yt)dt + a(y t )dB t <S> dy t = (^p{yt) - ^ '{yt)v{yt)\ dt + a(y t ) o dB t 

It immediately follows from equation ^ that any power of an iterated 
integral is a sum of iterated integrals and that a polynomial of iterated integrals 
is also a linear combination of single iterated integrals. For example: 

Y (l) Y (2,3) r (l,2,3) y (2,l,3) y (2,3,l) 

A O,T A 0,T ~~ ^0,T ^0,T ^0,T 

\^0,T J — ^0,T^0,T 

- oyM 
— 

^(0,1,0)^(1,1) _ V (0,1,0,1, 1) 9y (0,l, 1,0,1) ,0^(0,1, 1,1,0) 

A 0,T A 0,T ~~ A 0,T Z ^0,T J ^0,T 

^(1,0,1,0,1) 9y (1,0,1,1,0) ^(1,1,0,1,0) 
"r^O.T Z ^0,T ^0,T 



The number of terms from the product of iterated integrals grows fast, expo- 
nentially if all letters are different. 
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1.4 Expansions 

We are now in position to derive expansions for the solution of a differential 
equation. Picard iterations yield: 

y ,r(o,j) = 

Y 0>T (l,j) = [ f jl (Y s (0))dXP + ...+ [ f jn (Y s (0))dX^ 

J0,T J0,T 

= fji{Yo)X Q ^ T + ... + fj n (Y )X Q 7 l r 



(n) 



Y 0tT (2,j) = / f jl (Y s (l))dxW + ...+ f jn (y a (l))dXW 
Jo,t Jo,t 

= [ f jl (Yo, a (l)+Y )dXP + ...+ I f jn {Y s (l) + Y )dX ; 

J0,T J0,T 



Since io, s (l, •) is a sum of iterated integrals, the polynomials fji(Y s (l) + Yq) are 
also sums of iterated integrals and so is their integration with respect to X®. 
Yo,t(2) is thus a sum of iterated integrals and by recursion all Yo,t(?", •) are. A 
formal proof in the context of differential equations driven by rough paths can 
be found in [3]. 

Example 

Consider the Ornstein-Uhlenbeck process: dyt = a(l — yt)dt + bdWt and yo = 0. 
In this case, x[ 1] = t, xf } = W t and Y^ ] = for i 6 {1,2}. Applying Picard 
iterations, we obtain: 

Y , T (Q) = 



Y 0:T (l) = f o(l - 0)dXW + / 6dXj 2) 



r(l) ,7,^(2) 



r 0jT (2) = / T a(l - (aX« + bXP))dXP + /" MAf) 

= QjXq rjn d Xq rjp^ ClbX^ rp ^ ~\~ bX^ p 

^o,t(3) = clXq }p — o? ' Xq y ^ ~)~ T ' ^ cPXq ' ^ + h&Aq ^ ~\~ ^Aq ^ 

The solution of the stochastic differential equation can thus be approximated 
by a series of iterated integral of the drivers, whose coefficients are a function 
of the parameters. The iterated integrals capture the statistics of the drivers 
and are separated from the parameters. 

This derivation can be readily implemented in Mathematica ([6j^| (see [5] 
for an implementation of the shuffle product) but suffers a major drawback: 
each product of iterated integrals being a shuffle product, the number of terms 
produced grows extremely fast (exponentially in the worst cases) and rapidly 
becomes unmanageable. In the next section, we introduce a re-parametrisation 



■"■The code will work with Mathematica version 7, when parallel computing was introduced. 
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of the problem that circumvents this problem by providing an alternative rep- 
resentation of the expansion which can be processed in a distributed manner, 
alleviating large memory requirements. 

2 Re-parametrisation of the polynomials 

One thing to note in the derivation of expansions is that each successive itera- 
tions requires the explicit linear combination of iterated integrals for the previ- 
ous iteration; evaluating Yq^t{t) requires the complete expansion for Yqt(?— 1). 
As the number of terms grows, manipulating this object rapidly becomes un- 
wieldy. 

In this section, we describe a re-parametrisation of the polynomials which 
bypasses the need for an expansion in terms of iterated integrals. It provides 
a more compact representation of the approximate solution and is naturally 
amenable to parallel processing. For clarity of exposition, we first introduce 
the approach in the case of a one-dimensional (m = 1) differential equation 
with n drivers. We assume the polynomials fu to be of degree less or equal 
to q. In the last subsection, the procedure is generalized to m-dimensional 
differential equations. 

2.1 One-dimensional case 

We first remark that a polynomial P(y) can be written in terms of y — yo by 
writing its Taylor expansion around yo: 

p(y) = itli dkp (y°)(y-y°) k ( 9 ) 

fc=0 

Next we introduce the following new operation for iterated integrals: 

Kt > Kt = f KudKu = f xi u xl- u dxf (io) 

J s J s 

where /3- is the word (3 with the last letter removed and /3 en d the last letter 
of /3. This is a non-associative, non-commutative operation - in fact, it can be 
viewed as a non-commutative dendriform. We can rewrite the Picard iteration 
using this operation. Note that from now on, the interval [s, t] will be fixed to 



[0, T] and will be omitted. 



Y(0) = yo- yo = (11) 

Y(l) = [ T f 3 i(Ys(0) + y )dX^ (12) 
i=i ^ 

= ^/.i(yo)^ (i) (13) 
i=i 

n 

y(r + l) = ^2fii(Y a {r) + y ) > X« (14) 
i=i 

= jzjzh dkfMYs{r)k>x{l) (15) 

i=l fc=0 

= £ f((^«) fe > D^/iifc*)*®)! (16) 



fc=0 \ i=l 

Therefore, if we define the objects Q as: 



n 



i=i 



the Picard iteration takes the following form: 

y(i) = Q° 
<? 

y(r + i) = £y(r) fc >Q fc 

A:=0 

where y(r) fc is the usual product. 

2.2 Description of the approach 

Consider a one-dimensional differential equation with quadratic functions fu, 
i.e. q = 2. The new representation yields: 

y(i) = q° (is) 

y(2) = l>Q0 + Q0>Qi + (Q°) 2 >Q 2 (19) 

= Q° + Q° D> Q 1 + (Q ) 2 > Q 2 (20) 

y(3) = Q° + (Q° + Q° > Q 1 + (Q ) 2 > Q 2 ) > Q 1 (21) 

+(Q° + Q°>Q 1 + (Q°) 2 >Q 2 ) 2 >Q 2 (22) 

where (Q k ) e is the A:*' 1 object Q to the power Crucially, this new repre- 
sentation does not require the explicit computation of the shuffle products, 
keeping the number of terms under control. Moreover, the expression can be 
expanded into its summands and each summand be processed independently. 
Thus, we avoid the handling of a large expression and are able to parallelize 
the computation. 
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Note that the representation only depends on the maximum degree of the 
polynomials but not on the number of drivers. The Picard iteration can there- 
fore be done only once, stored in a file and used at a later date for any system 
that uses the same maximum degree q. 

Using this representation, it is possible to derive the stochastic expansions 
through a few stages: 

1. Expand the expression Y(r) into its mononials u. Each monomial u is 
a function of Q's that uses the non-commutative products. Importantly, 
the objects Q and the product > are only used as place-holders at that 
stage. For example, the first three monomials for Y(3) are Q , Q° > Q 1 
and (Q° > Q 1 ) > Q 1 . 

2. For each monomial u, the objects Q are replaced by their values from the 



model at hand. Each Q is a weighted sum of the drivers X® (eq 17), so 
each u becomes a polynomial V of in terms of the non-commutative 
product. As in the previous step, the product is still only employed as a 
place-holder and not instantiated. 

3. Each polynomial V is expanded into its monomials v. Each v is a function 
of X® and >. 

4. For each monomial v, the product is instantiated; its actual definition in 
terms of shuffle product is only used at this later stage. 

Each process only requires a fraction of the memory that a direct approach 
(replacing Q's by the X" and using the product's definition) would. Moreover, 
at all stages, each monomial can be processed independently from the rest, 
leading to natural parallelization. It is also important to note that the whole 
expression is actually stored in a file and thus not in memory at any point. The 
exact details of this procedure are given in section [3} 



2.3 Generalization 

In the multidimensional case, the Taylor expansion of a polynomial requires 
a larger number of terms that involve cross-products between the components 
of the vector Y. The objects Q are not indexed by k £ {0, . . . , q} but by 
the set OW m (0, q) of the ordered words of length up to q written with letters 
{1, . . . , m} and are now defined as: 

n III 

QrE^/ii^ (i) (23) 
i=i 

for j G {1, . . . , m}. The constant c(r) is the number of different words we can 
construct using the letters in r. The Picard iteration becomes: 

Y(l)O') = and Y(r + 1)^ = ^ Y(r) T > Q] 

r6OW m (0,<3) 



3 Implementation 



8 



This section describes how this new approach was implemented in Mathematica. 
In this part, iterated integrals of the drivers (X^ 1 >—>*"•)) are denoted j( J i>--->*™) 
to follow convention and drivers are numbered from with the first driver 
representing time. 



3.1 Shuffle product 

Since each product of two iterated integrals is a shuffle product, special care 
must be taken of its implementation. [5 J has shown a way of writing the product 
in Mathematica: 

Tocino [j [{x_}] , j[a_List]] : = 

Sum [j [Insert [a, x, k]], {k, 1, LengthOa + 1}]; 
Tocino [j [a_List] , j[{x_}]] := Tocino [j [{x}] , j[a]]; 
Tocino [j [a_List] , j[b_List]] := 

Ap [Tocino [j [a] , j[Drop[b, -1]]], LastOb] + 

Ap [Tocino [j [Drop [a, -1]], j[b]], Last@a] /; (LengthOa > 1 && Length@b > 1) ; 

This is a direct and natural translation of the following result: 
J a jP = [ J a -jPdJ a ™ A + / J a J^dJ^ d 



We present a new implementation of the shuffle product. The product is 
done iteratively instead of recursively and is based on string transformation. 
To calculate the shuffle product of two words a and /3, we first set the string 
'aa...aabb..bb', with as many a's and b's as there are letters in a and fi re- 
spectively (line 2). We then replace all occurences of 'ab' by 'ba' (line 6) and 
keep on iterating the replacement until none are possible (i.e when the word is 
'bb..bbaa..aa') while keeping track of the new words generated (lines 4-6). The 
set of all words thus created is the shuffle product of two arbitrary words of the 
same lengths of a and /3. Finally, the letters a and b in each word are replaced 
by their actual values from a and (3 (lines 7-9). 

1 Shuf f le [a_List ,b_List] : =Module [{list ,u, v} , 

2 u={StringJoin@Join[Table["a" , {Length [a] }] , Table ["b" , {Length [b] >]]} ; 

3 v=Table [0 , {StringLength [u [ [1] ] ] }] ; 

4 list=Flatten@NestWhileList [ 

5 DeleteDuplicates@Flatten@ 

e StringReplaceList [#, "ab"->"ba"]&,u,#!={}&] ; 

7 (v [ [# [ [1] ] &/@StringPosition [#, "a"] ] ] =a; 

s v [ [# [ [1] ] &/@StringPosition [# , "b"] ] ] =b ; 

9 v)&/@list] ; 



Figure 3.1 shows the relative performances of the implementation. The 
product of two iterated integrals with random words is calculated with the 
recursive definition and the iterative definition. The ratio of the time spent by 
the former over the latter is recorded and the process is re- iterated 1,000 times 
for each final-word length. On average, the iterative method is a bit slower on 
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Figure 1: Log2 of the average of the relative speed of the iterative over the recursive 
method. 

short words but faster for long words (from a total length of 8). Since longer 
words are more numerous, using the iterative method will be advantageous. 

3.2 Non-commutative product 

The non-commutative and non-associative product > is implemented as in 
Mathematica, an operation with no built-in meaning. We only set a few basic 
properties for this product: 

Unprotect [CircleDot] ; (* = esc c . esc *) 

10x_ := x; 

x_01 := 0; 

O0x_ := 0; 

x_0O := 0; 
Protect [CircleDot] ; 

Since is only used as a placeholder, its actual definition in terms of shuffle 



NCP[0, j [b_List]] := 0; 
NCP[1, j [b_List]] := j [b] ; 
NCP[j [a_List] , 0] := 0; 
NCP [j [a_List] , 1] := j [a] ; 
NCP [j [a_List] , j[{>]] := j [a] ; 
NCP [ j [a_List] , j [b_List] ] : = 



NCP [n_* j [a_List] , j [b_List] ] : = n*NCP [j [a] , j [b] ] ; 
NCP [j [a_List] , n_*j [b_List]] : = n*NCP [ j [a] , j[b]]; 
NCP[n_*j [a_List] , m_*j [b_List] ] : = n*m*NCP [ j [a] , j[b]]; 
NCP[x_ + y_, zj := NCP[x, z] + NCP[y, z] ; 
NCP[x_, y_ + zj := NCP[x, y] + NCP[x, z] ; 



product (eq. 




is coded in another function, NCP: 



Ap[j [a]*j [Drop[b, -1]], Last[b]] /; Length [b] > 0; 



10 



NCP[x_*(y_ + z_) , t_] := NCP[x*y, t] + NCP[x*z, t] ; 
NCP[(y_ + z_)*x_, t_] := NCP[x*y, t] + NCP[x*z, t] ; 



3.3 Picard iteration 

The usual Picard iteration is implemented with a helper function Picardlteration 
as following: 

Picardlteration [f _List ,X_] := 

Total [Maplndexed [ (Ap [#1 [X] * j [{}] , First [#2] -1] ) & , f ] ] 
Picard [f_List,XO_,n_Integer] := 

Nest [ (Picardlteration [f ,#])&, XO , n] 

For example, Picard [f ,x0, 4] outputs the stochastic approximation of the 
SDE with the functions / collected in a list in the first argument. This was 
used in [3] for a system with linear drift and quadratic variance. 

With the new representation in Q, it can be written directly as in eq. [2] 

PicardQlDim[Q_, R_, q_] := 

Nest[Q[0] + Sum[(#~r)0Q[r] , {r, 1, q}] &, Q [0] , R - 1] ; 

R is the number of iterations to be calculated and q the maximum degree of 
the polynomials /. PicardQIDim [] produces a very compact representation of 
the expansion, which needs to be processed further in order to give the same 
result as Picard [] . 



3.4 Distributed processing of monomials 

Going from the compact representation provided by PicardQIDim [] to the 
linear combination of iterated integrals j's is done in a few stages. Each stage 
modifies the representation of the stochastic expansion in such a way that a) 
computational requirements are minimized and b) it can be parallelized. 
As described in section |2T2| the workflow goes as follows: 

1. PicardQIDim [] produces a polynomial in Q and 0. 

2. Each monomial (in Q and 0) is extracted and stored in a list (actually a 
file). 



For each monomial, the Q's are replaced by their values in j (eq 17). They 
are now polynomials in j and 0. 

The monomials (in j and 0) from each polynomial are extracted and 
stored in a file. 



5. The product is instanciated in terms of the shuffle product (eq. 10). 
The result is a linear combination of iterated integrals j for each mono- 
mials. The stochastic expansion is the overall sum of the all those linear 
combinations. 



As can be seen on figure 3.4, the different polynomials are successively 
expanded in terms of monomials, which in turn are processed independently. 
Since each expression is effectively broken down in small parts and dumped 
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Stage 1 


(Q° + Q 1 &Q 1 )QQ 2 ... 





Stage 2 


((2Q° Q Q 1 ) Q Q 1 ) Q Q° 









Stage 3 (ajW + . . .) (a 2 «/« + . . .) 



Stage 4 (( a jW) (a 2 J«)) 



Stage 5 



_J_ 3 & 5 7(1,2,1,1) 4. 
4080 " J ^ 



l a 2 6 5j(2,2,l,l) 



Figure 2: Workflow. From the compact representation in Q and to the linear 
combination of iterated integrals J's. 

into a file, a much larger of terms can be computed, as can be seen in the next 
section. 

Since each processing stage follows roughly the same logic, with slightly 
different transformation rules, we only detail the first one: going from the 
polynomial in Q and to its monomials. This is done step by step, and not by 
applying a global rule, since doing so would produce a very large output and 
run the risk of failling due to memory limits. 

1 Print ["=> Deriving the monomials in Q and 0"]; 

2 changeHappened = True; 

3 workingFile = "uexpansion" ; 

4 finalFile = "uexpansion_Final" ; 

5 DistributeDef initions [workingFile , finalFile, rulePlus, rulePower ,CircleDot] ; 

6 (* write expansion in *) 

7 in = OpenWriteSmainFile ; Write[in, PicardQlDim[Q, 4, 2]]; Close@in; 
s SetSharedVariable [changeHappened] ; 

9 ParallelEvaluate [ Put [finalFile <> ToString [$KernelID] ] ; ] ; 

10 While [changeHappened, 

11 changeHappened = False; 

12 positionList = GetStreamPositions@workingFile ; 

13 DistributeDef initions [positionList] ; 

14 ParallelEvaluate [in = OpenReadOworkingFile] ; 

is ParallelEvaluate [Put ["temporary_" <> ToString [$KernelID] ]] ; 

16 WaitAllOTable [ 

17 With[{sp = sp}, 
is ParallelSubmit [ 

19 SetStreamPosition [in, sp] ; 

20 v = Read [in, Expression] ; 
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21 If[FreeQ[v, Plus] && FreeQ[v, Power], 

22 PutAppend[v, finalFile <> ToString [$KernelID] ] , 

23 If [Head [v] === Plus, 

24 changeHappened = True; 

25 Scan [Put Append [#, "test_temporary_" <> ToString [$KernelID] ] &, v] , 

26 firstPlus = MinSMap [Length®* &, Position[v, Plus]]; 

27 firstPower = Min@Map [Length@# &, Position[v, Power]]; 

28 If [{firstPlus , firstPower} != {Infinity, Infinity}, 

29 If [firstPlus <= firstPower, 

30 changeHappened = True; 

31 Do [v = Replace[v, rulePlus, depth], {depth, firstPlus - 2, 1, -1}] ; 

32 v = Replace [v, rulePlus]; 

33 Scan [Put Append [# , "temporary." <> ToString [$KernelID] ] &, v] , 

34 changeHappened = True; 

35 v = Replace [v, rulePower, {firstPower - 1}]; 

36 PutAppend[v, "temporary." <> ToString [$KernelID] ]] , 

37 PutAppend[v, "temporary." <> ToString [$KernelID] ]]]]]] , {sp, 

38 positionList}] ; 

39 ParallelEvaluate [CloseOin] ; 

40 ConcatenateFiles [tempFileNames , workingFile] ] ; 



4i ConcatenateFiles [Map [finalFile <> ToString@# &, kernellDList] , finalFile]; 

The data to process is in workingFile. At each step, an entry is read from 
this file and processed by one of the kernels (lines 20-37). If the result is a 
monomial, it is added to the final file (line 22). Otherwise, it still requires 
further processing and is added to a temporary file (lines 25,33,36,37). When 
all entries in workingFile have been processed, workingFile is replaced by 
the temporary file (line 40). The algorithm loops back until all monomials have 
been extracted to finalFile. Note that each kernel has its own temporary and 
final files, as writing to the same file in parallel typically results in corrupted files 
and missing entries. Processing an entry consists in finding the highest 'Plus' 
in the expression and distributing the monomials around it. In this manner, 
the polynomial is first decomposed into the largest polynomials, which are 
processed further separately. 



4 Expectation of an iterated integral 

It is often of interest to calculate the moments of the solution of the SDE 
and this can be approximated by computing the moments of the stochastic 
expansion. Since the expansion is a weighted sum of iterated integrals j, its 
expectation is simply the weighted sum of the integrals' expectations. If the 
drivers consist of time and Brownian motions, the expectation of an iterated 
integral has a simple analytic form that can be arrived at recursively ([5]). 

Here we present a more direct way of calculating this quantity. Given a 
word a and assuming Wiener processes ([2]): 

Eja(t) / ^ ^ a 15 a se 1 uence °f an d pairs mm 
I Pa^—\ otherwise 



where p a =2 2 and q a = n(#{"i / 0}) + = 0}). 



2 ""^ — 2 

Thus, for example: 

£j(o,i,i,o,o) ' = i/2 2 /¥ 3 + 2 / 2 )/(3 + 2/2)! = g 

^(0,1,1,0,0,1) = o 

£j(2,2,l,l,3,3) = l/26/2 t (0+6/2) /(0 + 6/2)! = J 

£j(2,2,0,l,l,3,3,0,0,0) = l/ 2 6/2i(4+6/2) /(4 + 6/2) , = 

This result is implemented in Mathematica. The expectation for a word a 
is then calculated in at most \a\ steps: 

ExpSBM[t_, j[a_List]] := Module [{i, c}, 
i = LengthSa; 
c = {0, 0>; 
Catch [ 
While [i > 0, 
If[a[[i]] == 0, 
c += {0, 1}; i— , 

If[(i > 1) && (a[[i]] == a[[i - 1]]), 
c += {1, 1}; i -= 2, 
c = {Infinity, 0}; ThrowOO 
]]]]; 

(1/2) ~First@c t~Last@c/(Last@c) !] ; 



5 Example 

Consider the following SDE: dY t = a(l - Y t )dX x + bY t 2 dX 2 and Y = 0. In this 
case, m = 1, n = 2, q = 2. The two functions / are fi,i(x) = a(l — x) and 
/i,2 0) = bx 2 . 

Only two things are required from the user: the definition of the objects Q 
in a transformation rule, easily obtained by derivation (eqfT7|), and the number 
of Picard iterations to be computed. In this case, the three Q's are: 

Q° = ^(a(l-0)X( 1 + b0 2 X^) 

= aXW 
Q 1 = i(-aX( 1 ) + 260X( 2 )) 

= -aXW 
Q 2 = i(0X( 1 + 26X( 2 )) 
= bX^ 



Therefore, the transformation rule corresponding to this system is: 

ruleModel={Q [0] ->aj [{1}] , Q [1] ->-aj [{1}] , Q [2] ->b j [{2}] }} ; 

If the first driver is time, we can follow the convention that drivers are numbered 
from 0: 
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ruleModel={Q [0] ->aj [{0}] , Q [1] ->-aj [{0}] , Q [2] ->b j [{1}] }} ; 

A small number (e.g. 4) of Picard iterations is sufficient for a good approxi- 
mation (see [1]). An optional parameter is the maximum word length of the 
iterated integrals. 

Running the algorithm with these parameters, we obtain the expansion, 
which is a weighted sum of 676 iterated integrals and already has a ByteCount 
(size) of 394'288. Its expectation is a much smaller expression: 

aT - a -T 2 + — T 3 + (-a 3 b 2 - -)T 4 - -a 4 b 2 T 5 + — a 5 6 2 T 6 + (— a 5 6 4 - -a% 2 )T 7 
2 6 4 24 ; 20 360 v 140 24 ' 

. 1 7 ,2 21 fif4Nm8 157a 7 6 4 T 9 A3a 7 b G 17a 8 5 4 in 1 „, fi ,, 

+( a 7 b 2 a 6 6 4 )T 8 + + )T 10 a 8 6 6 T n 

v 192 160 ; 3024 v 1800 2800 ; 100 



This is confirmed by the previous implementation of Picard iteration with 
the simple code: 

P0[x_] := a (1 - x) ; 
Pl[x_] := b x~2; 

expansion = Picard [{P0, PI}, 0, 4] ; 
ExpSBM [t , expans i on] 

However, if we now set the initial value to an arbitrary yo instead of 0, the 
stochastic expansion contains a much larger number of terms. It is calculated 
in the same manner using Q as: 

ruleModel={Q[0]->a(l-yO)j [{0}]+by0~2j [{1}] , Q[l]->-aj [{0}]+2by0j [{1}] ,Q[2]->bj [{1}] }; 

This results in an expansion with 10'710 unique iterated integrals and a ByteCount 
of 15'841'624. This cannot be confirmed by the simple code as it runs out of 
memory and crashes before completing. Note that the expansion is usually left 
in a file whose entries are parts of the linear combination. Thus, it is possible 
to, for example, compute the expectation of the expansion without having to 
store it in memory at any point. Moreover, while computationaly expensive, 
these expansions can be calculated once in a general case and saved in a file for 
future application without having to recompute them de novo. 



6 Conclusion 

Stochastic expansions provide a local approximation of the solution of a stochas- 
tic (or deterministic) differential equation. They can be used for a variety of 
applications, from simulation to parameter estimation. However, as the num- 
ber of terms grows exponentially with the desired precision, they can rapidly 
become unwiedly to manipulate. 

We presented a new way of calculating these expansions that bypasses the 
limitation of the usual approach, via a re-parametrisation of the problem and 
the parallelization of the computation. We have shown that in a simple example 
our method was able to compute the expansion when a direct approach failed. 
We also presented two new approaches for efficiently deriving the shuffle product 
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of two iterated integrals and the expectation of an iterated integral, when the 
drivers are time and Brownian motion. 

So far, our approach has been implemented for one-dimensional differen- 
tial equatiorj^j However, the theoretical foundation for the multi-dimensional 
case is available, as presented in section 2.3 Now that the computing require- 
ments have been alleviated, an implementation for the general case is possible. 
Stochastic expansions will then be available for more complex systems. 
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